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Abstract 

In this paper we calculate the ground state phase diagram of a Josephson 
Junction ladder when screening field effects are taken into account. We study 
the ground state configuration as a function of the external field, the penetra- 
tion depth and the anisotropy of the ladder, using different approximations 
to the calculation of the induced fields. A series of tongues, characterized by 
the vortex density u>, is obtained. The vortex density of the ground state, 
as a function of the external field, is a Devil's staircase, with a plateau for 
every rational value of u. The width of each of these steps depends strongly 
on the approximation made when calculating the inductance effect: if the 
self-inductance matrix is considered, the uj = phase tends to occupy all the 
diagram as the penetration depth decreases. If, instead, the whole inductance 
matrix is considered, the width of any step tends to a non-zero value in the 
limit of very low penetration depth. We have also analyzed the stability of 
some simple met ast able phases: screening fields are shown to enlarge their 
stability range. 



PACS numbers: 74.50.+r, 05.20.-y, 64.70.Rh 



I. INTRODUCTION 

Theoretical research in Josephson junction arrays (JJA) is continuously progressing 
through models which involve an increasing complexity. They represent a better approx- 
imation to the understanding^ and prediction of the many different interesting phenomena 
which occur in such systems.^ An important contribution to this advance is the inclusion 
of current induced magnetic fields (CIMF) developed by different groups in the last years.ll 
Taking CIMF into account is compulsory in order to provide a correct description of Joseph- 
son junction arrays at low temperature, when the penetration depth of the magnetic field is 
of about the cell size. 
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In all the cases, the study was carried out through the numerical simulation of the 
dynamics of the gauge invariant phase differences. Interest has been mainly focused on the 
effects of CIMF in the properties of arrays driven by external currents. However, we have 
no knowledge of studies on the ground state properties of inductive arrays. 

This paper deals with the static properties of a Josephson junction ladder (JJL), with 
anisotropy in the Josephson couplings, in the presence of an external magnetic field (figure 
|TJ). In particular, we have calculated the ground state phase diagram of a system defined 
by a hamiltonian which includes the magnetic energy due to the CIMF, in addition to the 
usual Josephson coupling contribution. 

Recently two different groupsHi have faced the ground state problem of a JJL in the 
presence of a magnetic field and in the limit of infinite penetration depth - no screening 
effects -. In this approximation, some general properties concerning the ground states can 
be deduced. In regard to the ground state problem, the hamiltonian describing the system 
belongs to a universality class of convex 1-D models of spatially modulated structuresjl 
such as the Frenkel-Kontorova model.! This fact allows an exact description of the ground 
states phase diagram and the relevant elementary excitations of the systemJl The diagram, 
a function of the external field and the anisotropy parameter, consists of a series of tongues 
labelled by the vortex density uj. Both rational and irrational values of to are possible, 
corresponding, respectively, to commensurate and incommensurate phase configurations. 
The vortex density of ground state configurations as a function of the external magnetic field 
is a Devil's staircase function, with plateaus for every rational value of to. Incommensurate 
ground states exhibit two regimes, separated by an Aubry transition!: below a certain value 
of the parameter that describes the anisotropy, the configuration is undefectible (no defects 
can be sustained) and unpinned (any external current, though infinitesimal, causes a non- 
zero voltage); above this value, the solution is defectible and pinned. 

In this paper we use the work by Mazo et al.l as starting point and include CIMF. We thus 
obtain a more realistic description of the ground states and, in general, of the equilibrium 
properties of the ladder. Such results may be of interest in understanding experiments in 
JJL where relevant parameters can be fixed at will. 

We have used different numerical methods - effective potentials method combined with 
root finding methods and stability analysis of solutions, as well as dynamical relaxation - in 
order to find the ground state phase diagram and other metastable configurations. Results 
are obtained for three approximations to the calculation of the induced fields: (A) self- 
inductance contributions; (B) self-inductance plus nearest neighbours mutual inductance 
contributions; (C) full-inductance matrix. In all the cases the vortex density of the ground 
state, as a function of the external field, exhibits a Devil's staircase structure. Special 
attention has been focused on the behaviour of the system in the small penetration depth 
limit. In this limit, a vortex can be described as a flux quantum concentrated in one only 
cell. The ground state phase diagram shows an important difference depending on the 
approximation made when calculating the inductance effect: if the self-inductance matrix is 
considered, the uj = phase, in this limit, nearly occupies all the diagram. When the whole 
inductance matrix is considered we find ground states with no null vortex density in a wide 
region of the phase diagram. 

The variation of the induced flux with the penetration depth allows an estimation of 
the physically interesting range of values of this array parameter. We have also studied the 
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dependence of some of the vortex properties with the penetration depth and the anisotropy 
of the ladder, such as its extension and the distribution of gauge-invariant phases and the 
induced flux. Finally, we consider the stability of some simple metastable commensurate 
phases when the external field is varied. Notably, the stability intervals enlarge when the 
penetration depth decreases, existing a critical value of the penetration depth for the stability 
of each phase at zero external field. 

The paper is organized as follows: in section II we introduce the model and the different 
methods and approximations used to compute its properties. Results on the ground state 
phase diagram and the stability of some simple commensurate phases are reported in section 
III. Different approximations to the calculation of the CIMF are discussed. 



II. DESCRIPTION OF THE MODEL AND THE METHOD 

The classical hamiltonianB describing the system is: 

H = - Ei [J x cos(9i - 9 i+1 - tt/o - tt/<) 

+ J x CO S (6' i - 9' i+1 + 7T/o + Ttfi) + JyCOS^i - 9'i)] (1) 
+ 2^0 fiLij fj- 

Here 0j, (6^) denotes the phase of the superconducting order parameter on the upper (lower) 
branch of the ladder at the ith site (see fig. [I]), /q is the magnetic flux due to the external 
field, which is assumed to be constant along the array, fi is the induced flux through 
plaquette i, a function of the currents in the ladder. Both /o and fi are expressed in terms 
of the flux quantum, $ - Thus, the total magnetic flux through a given plaquette i 
is = $ ex * + $* nd = $q(/o + fi)- The model is periodic in f with period 1 and has 
symmetric reflection about f = | in the interval [0, 1]. Thus we will restrict our analysis to 
values of f within the interval [0, |]. 

When writing Eq.([l|) we have made a convenient gauge choice: we consider that the 
vector potential is parallel to the long axis of the ladder, and takes opposite values on the 
upper and lower branches (see Fig.[l|). In this gauge fo and fi are trivially related to the line 
integrals of the vector potential a through a link of the ladder: a a p — f£ adl = e7r(/ + /»), 
where e = +1(— 1) for upper (lower) links in the ladder and e = for vertical links. 

J a [a = x,y), the Josephson coupling energy, is related to the critical current through 
the junction, I ca , by J a = J ca $ /(27r). The inductance matrix of the array, L, is defined as 

where Ay is an adimensional matrix containing just geometrical coefficients (see Appendix 
A). A_i_ is the penetration depth, defined as in! 

> - 1 *° (3) 



47r 2 fi J x a 
being a is the lattice spacing. 
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It can be seen that the configurations which minimize hamiltonian ([!]) comply with 
9i + 9 i = const. Fixing this constant equal to and normalizing by J x in order to work with 
adimensional quantities we get 

H = - Ei [2cos{Bi - 9 l+1 - nf - tt/0 

where the quotient J y J J x defines the anisotropy of the ladder. To solve the ground state 
phase diagram we will restrict our analysis to expression (|J) . 

We consider here three approximations to the inductance matrix: the simplest model 
(case A) assumes a diagonal inductance matrix. In this case, the flux induced in a given 
plaquette only depends on the mesh current in the same plaquette. Next step in complexity, 
case B includes also nearest neighbours inductances for the couplings between cells. Then, 
we assume L^ 1 = L5ij + M8ij±\. Case C considers the full-range inductance matrix. In the 
first case the term in Eqs.(|l]) and giving an account of the magnetic energy becomes 

Hmagn = J2d K fi, with d K = 87T 3 A_l(A _1 ) 

i 

In case B, 

Hmagn = d xfi + ad Kfi(fi-l + fi+l)- (6) 

i i 

In the system under simulation (we consider square cells; currents are supposed to flow 
within a cylinder of length a and ratio 0.005a) dx ~ 6.8Aj_ and a = %r ~ 0.21355. 

We have used different methods to solve the problem. In the cases A and B it is possible 
to do it using an effective potentials method properly adapted to study our model,El which 
is numerically equivalent to a ID system with just next-nearest-neighbours interactions. 
Effective potentials method is an efficient method to study the ground state configurations of 
such kind of systems in the thermodynamic limit. This method is based on the computation 
of certain functions, the effective potentials, which contain all the relevant information on 
the relaxation of local fluctuations to the ground state configurations. A long computation 
time is required if one wants to obtain the phase diagram of the system with a high precision. 
This suggests the convenience of complementing the method with other procedures. 

Effective potentials can provide, within a reasonable amount of time, approximate solu- 
tions to the ground state problem as a function of the external field, the anisotropy of the 
ladder and the penetration depth. Starting from these guesses, one can obtain more precise 
results by applying standard root finding methods (calculating stable solutions to |^ = 0) or 
even dynamically (letting the approximate solution relax). We make note that root finding 
methods require the use of Eq. ([[]) to describe the system since Eq. (|) is just an adequate 
expression when dealing with minimum energy configurations, which is a reduced subspace 
of the whole system. We have checked that the same results are obtained if one applies 
the effective potentials method with a high precision or if one combines it with any of the 
complementary procedures described above. By comparing the energy curves corresponding 
to different configurations one can determine the border of the tongues with different vortex 
densities. 
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(4) 



(5) 



Moreover, making use of these procedures one can study how a ground state configuration 
modifies when the parameters vary. In this case, it is convenient to keep in mind that, in 
general, a vortex configuration is stable beyond the range of parameters in which it is 
the ground state solution. There, root finding methods are adequate, and they must be 
completed by doing the linear stability analysis of the solutions. 

Model C involves the total flux matrix. Interactions between the variables extend to 
all the lattice and the problem can not be tackled with the effective potentials method. 
In this case, we consider the results achieved in approximation B as guesses and let the 
system evolve dynamically and relax down to the equilibrium. Details about the dynamical 



Figure §a shows the ground state phase diagram in the case of infinite penetration depthEl 
Aj_ (thus neglecting CIMF). The different tongues are characterized by the vortex density uj. 
This quantity is directly related to the periodicity of the configuration: a value uj = | implies 
that relevant physical quantities - gauge invariant phases differences, induced fluxes... - are 
spatially periodic: every q plaquettes these quantities are exactly repeated. Here vortices are 
defined as usual. We make use of the well-known property of fluxoid quantization to define 
the vorticity n p on any plaquette. The clockwise sum of the gauge invariant phases (restricted 
to the interval (— ir, ir]) along the links of the cell gives Z) Q /3gi(@a — 0/3 — ^a/?) = 27r(n p — f tot ). 
The vortex density u is equal to the spatial average of n p . 

As mentioned in the Introduction, the vortex density as a function of the external field 
is a Devil's staircase, with plateaus for every rational value of uj. Figures 0b and ^c show 
the phase diagram in the cases A and B, for a penetration depth Aj_ = 1, computed using 
the effective potentials method. Diagrams (b) and (c) are qualitatively similar to diagram 
(a). As one expects, the CIMF tend to push the external magnetic field out of the array: as 
Aj_ decreases, the u = tongue grows. There is, however, a remarkable difference between 
these diagrams: while in case A all the tongues but the uj = one compress, in case B the 
uj — | phase does not shrink, and the rest of the phases are compressed between these two. 
In short we will study carefully the limit of this behaviour when Aj_ —>■ 0. 

The effect of the CIMF is to increase the critical field f c up to which the uj = configu- 
ration is the ground state. The devil's staircase is thus restricted to a narrower interval of 
values of the field. 

The question arises whether CIMF are able to change qualitatively the nature of the 
phase diagram. In other words, we want to check if for some Aj_ the array is able to expulse 
completely the external field and, consequently, this tongue occupies all the phase diagram. 
If this is not the case, is the devil's staircase structure preserved for all the values of Aj_? 

In order to gain a more complete understanding of the properties of the model and, 
in particular, to throw light on the previously raised questions, it is interesting to study 
the dependence of the vortex characteristics on the different physical parameters. Such 
study is firstly carried out by considering commensurate phases with a low vortex density 
(e.g. uj = Y^g) in order to prevent vortex-vortex interaction effects. In particular, we are 
interested in studying the vortex extension. It is directly related to the distribution of the 



algorithm are given elsewhere. 



III. RESULTS: GROUND STATE PHASE DIAGRAM AND STABILITY 

ANALYSIS 
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gauge invariant phases around the vortex barycenter and depends on the values of J y and 
\± (see figure §). in cells near the vortex centre, the phase decays exponentially. This decay, 
as distance increases, becomes smoother; for large % (where i is the distance to the centre), 
the phase is of the form ~ i~ 3 ). Anisotropy affects the exponential part of the curve: 
a decrease of J y implies a smoother exponential decay, while the long-distance behaviour 
remains unchanged. Instead, varying Aj_ makes the whole curve shift. In cells far away 
enough from the centre, the flux is negative and its absolute value is a decreasing function 
of the distance. The negative flux is due to the sign of the mutual inductance term (see 
Appendix A). This feature was already reported by Phillips et al.i. As Aj_ decreases the 
vortex becomes more and more localized on its central plaquette and, in the Aj_ — > limit, 
tends to identify with the fluxoid (the quantized magnitude defined above, that indicates 
the barycenter of the vortex). We can think of Aj_ as the radius of the vortex: Aj_ — > 
implies that the vortex remains restricted just to the cell where n p = 1. 

When no CIMF are considered (Aj_ — > oo limit) the total magnetic flux through a plaque- 
tte is just the external flux, which is constant along all the array. Thus, the flux distribution 
along the ladder is independent on the vorticity. Then, there is no flux quantization and the 
vortex density is not a flux quanta density but a fluxoid quanta density. On the contrary, 
the situation changes drastically when CIMF are taken into account, being the number and 
extension of vortices directly connected to the distribution of the induced flux along the lad- 
der. Let us consider the limit Aj_ — > and the lj = ground state configuration; there, the 
currents tend to uniformly screen the external field (in every cell fi — > — /o). Thus, the array 
exhibits a behaviour that resembles the Meissner effect: the external field is screened by the 
system and no field penetrates the array. For any other value of u, the flux distribution is 
quite different. In the cells where the vorticity is equal to zero, the induced flux tends to 
cancel the external one, being the total flux equal to zero. In the cells where the vorticity is 
equal to one, the induced flux fi — > (1 — / ), being the total flux equal to one flux quantum. 
Then, in the Aj_ — > limit the fluxoid identifies with fluxon, and it is well localized in a cell 
of the array. This is illustrated in figure |4], where the dependence of the induced flux on Aj_ 
is shown both in cells with vorticity and 1. We have chosen a configuration with oj = |, 
and fo — h, but the behaviour is general: for other values of u the fluxes in cells with zero 
vorticity depend on the distance to the nearest vortex, but this difference is of the order of A^_ 
in the Aj_ — > limit. We can distinguish 3 regions: For A^ > 4, \fi\ < 0.1 f , and considering 
an infinite penetration depth is a justified approximation. On the other hand, for Aj_ < 0.12, 
we observe the low Aj_ behaviour: > 0.9 (n* — fo), and screening fields are dominant. 
Between them there is an intermediate region, around A^ = 0.7 (where the derivative of the 
induced flux respect to the logarithm of the penetration depth is maximum). 

The description presented in the previous paragraph permits to obtain some simple 
expressions for the energies of the different configurations for low values of Aj_. The hamil- 
tonian ([]]) consists of two components, corresponding to the Josephson and the magnetic 
energies. We have numerically checked that, as Aj_ — > 0, the first term saturates first 
than the magnetic one. Thus, for low enough Aj_, we can approximate the energy per pla- 
quette by E p = -3 + (rij - / + 6 ft) (L^ 1 )^ (rij - f + Sfj), where 5 ft ~ O(Ajl), 

N p is the total number of cells, and n, = 0, 1 is the vorticity of cell i. Let us con- 
sider first the case A in this approximation. The energy per cell of a configuration with 
u = f is E p = -3 + d K (| (1 - f ) 2 + fo 2 )- As f E [0, |], (1 - fo) > fo, and 
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Ep is a increasing function of uj: for any fo ^ \, oj q < UJ\ implies E(uj ) < E(uj\) 
and the configuration with uj = is the ground state. If f = ~, as previously 
defined, has the same value for all u/s. A second order approximation in dx is re- 
quired. It is easily obtained that E{uj = ~; fo = ~) < £7(a; = 0; / — \) (i n particular, 
£(w = 0;/o = |)- E(u = §; / = |) = ^/(6tt 2 J,)). 

We remark that these results correspond to the case A and an extreme limit (Aj_ — > 
0). Nevertheless, the devil's staircase is observable down to low values of Aj_ (we have 
checked this point even at Aj_ = 0.012). As an example, figure [5] shows the energy of stable 
configurations with different values of uj (0, |, |, |, |, |) as a function of the external field 
when Xj_ = 0.5. The energy of the ground state staircase corresponds to the envelope of the 
curves, thus an approximation to the uj(fo) function can be obtained from them. 

Things change when one considers a more complete approximation to the inductance 
matrix (models B and C). In this case the uj = phase does not fill the diagram at any value 
of Aj_ and other commensurate phases are clearly appreciated. Before performing a rigorous 
analysis, we present a plausibility argument in support of this statement. Let's begin by 
comparing uj = and uj — ^ phases. In a uj = configuration, the currents and flux are 
identical in all the cells. On the other hand, a configuration with uj = | exhibits a spatial 
periodicity with period 2a; when f = | the flux and currents on one cell are of the same 
module and the opposed sign respect to those on the adjacent cells. This allows us to define 
an effective d K for both configurations, so that the magnetic energy per cell is just dxeffff 
(see Appendix A). In case B dxeffi^ = 0) = 10,945Aj_ and dxeffi^ — §) = 4, 617A^ while 
in case C dxeffi^ = 0) = 11, 176Aj_ and dxeffi^ = f) = 4, 638Aj_. In general, for wo < W\ 
dK e ff( w o) > d-Keffiwi), and thus E(wo, fo = 1/2) > E(wi, fo = 1/2). As the energy is 
a continuous function of f the previous inequality maintains for a range of / values near 
fo = 1/2. 

Moreover, it is possible to extend in a trivial way the argument previously developed for 
case A, and to calculate the energy per plaquette of any vortex distribution as a function of 
fo. Note that a configuration with uj = 2 is described by a periodic spatial structure the basic 
sub-array of which consists of q plaquettes containing p vortices. We can thus reduce the 
system to one with just q cells. In order to do that it is necessary to generalize the previous 
reckoning and to redefine the components of the matrix L in order to take into account 
the contribution of each of the infinite replicas of the basic sub-array. Thus the inductance 
between two cells at distance j is given by Lj = J2 n Lo,j+nq, where n = 0,±1,±2, 
j = 0, . . . , q— 1 and Lo,i — A),-i- 111 the limit Aj_ — ► the induced flux in any cell is given by 
a vector F = {f a fbfc---} with = (1 — f + 5fi) or (—fo + 8 fi) depending on the occupation 
number of the cell, rij. Thus, the energy per plaquette, up to the first order in Aj_, is given 
by E p = — 3 + 2J ° N fi(L )ijfj, with /j = — fo. We have computed this expression for 
a series of values of uj. By comparing the energies of the curves for different configurations 
we have obtained a Devil's staircase, see figure |^. This figure corresponds to the case C; in 
case B an analogous behaviour is observed. 

Hereafter, we will consider the response of the JJL to continuous variations of the exter- 
nal field. We will restrict our analysis to the study of the stability intervals of some simple 
commensurate phases which are the ground state solution at some value of the parameters 
of the system (thus, we consider only ordered phases including just vortices (and no antivor- 
tices) and for which < uj < 5). Such perspective, in the no screening field effects limit, has 
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been studied in reference! in order to characterize the dynamical approach to equilibrium, 
which has been shown to lead to slow relaxation. Here we will just focus on an important 
difference which appears when CIMF are considered. Such study is carried out through a 
quasi-static computation of ordered stable configurations (local minima of hamiltonian (fj)) 
with a determined vortex density, when the external field is slowly varied. 

As previously mentioned, the range of stability of any vortex configuration (u) is broader 
than the interval of the parameters in which it is the ground state. In general, there exists 
a critical value of the field for the stability of each phase: for f < f c {uj) the phase uj is 
no more stable. The loss of stability when f decreases occurs in this way: decrease in the 
external field implies an increase in the supercurrent through the horizontal links in order to 
maintain the vortex density in the array. The instability of the state takes place when the 
supercurrent in one link reaches its maximum value. At this point any small change in the 
field can not be sustained by an increase of the currents and the vortex structure becomes 
unstable, and the system relaxes to a new vortex configuration. That would be the process 
for the changes of vortex density when the external field is varied or when the thermal noise 
is high enough to produce a vortex to jump over the energy barrier of the metastable phase 
and then the system approaches to some other more stable phase. 

In the limit of neglecting screening effects! f c (uj) > for all uj 7^ 0, and thus when 
/o = only the uj = phase is stable. The inclusion of CIMF changes this situation. As Aj_ 
decreases the range of stability of a given phase enlarges. Moreover, for each configuration w 
there is a critical value for the penetration depth Aj_ c (u>): if Aj_ < Aj_ c (w) the phase is stable 
at every value of the external field. Let's consider the two extreme cases: the configuration 
containing one single vortex and the w = 1/2 phase. The repulsive character of the vortex- 
vortex interactions implies that the stability of the configuration containing a single vortex 
is a necessary condition for the stability of any phase with < uj < |, so that the particular 
value of A^ at which stability of the configuration occurs (A^_(/o)) is an upper bound for the 
stability of the phases with < uj < |. On the other hand, stability of the uj = | phase 
ensures the stability of any other phase with < uj < |, so that the particular value of Aj_ at 

which stability of the uj — | occurs (A^ 2 (/o)) is an lower bound for the stability of the phases 

with < uj < \. Thus, AY c 2 (/o) < X± c (fo) < Kdfo)- Figure ^ shows the regions of stability 
of the vortex configurations. For values of the parameters above the curves (region S), any 
vortex configuration is stable. In region S-I, as we move towards the origin of coordinates, 
the different states become unstable (in the order of decreasing w). In I the only stable 
configuration is that with w — 0. Looking again at the supercurrents in the array, we see 
that as Aj_ decreases the gauge invariant phase differences of the metastable configurations 
approach to zero and thus the supercurrents are lower, rendering the phase more stable. 



IV. DISCUSSION 

In section III we have presented the phase diagram of a Josephson junction ladder in 
the presence of screening magnetic field effects. Numerical results evidence the existence 
of a series of tongues labelled by the mean vorticity uj. Such magnitude exhibits a devil's 
staircase structure when the external field is varied. 

When CIMF are considered, the system presents a behaviour that resembles the Meissner 
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effect: the self-induced field tends to push the external field out of the array, causing the 
growing of the u> = tongue and the shrinking of the range of parameters where the devil's 
staircase is observed. We have compared the results of using different approximations when 
calculating the induced fluxes. If only the self- inductance term is considered, for a value 
of Aj_ low enough the ui = phase occupies all the phase diagram except for a tiny region 
near the fo = \ line. However states with inserted fluxons are stable and, moreover, their 
range of stability increases as Aj_ decreases. Instead, if one takes into account the whole 
inductance matrix, the critical field for the u = phase (the frustration above which the 
configuration is no more the ground state) remains lower than | for all the values of the 
penetration depth. Thus commensurate phases with vortices are always clearly visible in 
the phase diagram. 

The three approximations made: A (self-inductance), B (self- inductance plus nearest 
neighbours mutual inductance terms) and C (full-inductance matrix) correspond to different 
distributions of the relative weights of the inductance matrix components. We remark that 
these distributions are a function of the geometry of the currents flowing in the array (see 
Appendix). Let's consider the case in which currents flow inside cylindric tubes of radius 
r and length equal to the lattice spacing a (the qualitative conclusions can be extended to 
any kind of cross section). If r << a, > |A i)i+1 | >> \A iti+ j\(j > 1) and approximation B 
(considering just the self-inductance plus the nearest neighbours terms in A) is justified. As 
r increases, the terms A i>i+ j(j > 1) also increase, but are yet too small. They give just small 
corrections to the final results. In a narrow range of r values around r ~ 0.25 the dominant 
contribution to the inductance is self-term (and A is a good — th order approximation). 
Finally, for greater r, lA^+jl/A^ cannot be neglected and considering the whole inductance 
matrix is compulsory. 

This behaviour (the existence of an infinite set of ground states as the parameters vary 
which show a Devil's staircase structure) is characteristic of a broad class of spatially mod- 
ulated structures with convex interparticle interactions. In the limit of neglecting screening 
field effects (Aj_ — > oo) it is well established the equivalence, regarding the ground state 
problem, of hamiltonian (P with a one-dimensional XY model with anisotropy and the 
ground state problem of the system is equivalent to the one of a Frenkel-Kontorova model 
with convex interparticle interaction, which allows for applying the Aubry theory for this 
class of models. 

However, and despite the qualitatively similar behaviour shown by our simulations, the 
inclusion of CIMF renders it difficult to establish any equivalence between model (|I|) and 
the 1-D models named above. This point is beyond the scope of this paper and remains as 
an open question deserving future research. 

The introduction of the CMIF allows to study the continuous variation of the system 
from the full penetration of the external field (Aj_ — > oo; = f , for Aj_ > 4) limit to the 
case of an array where all magnetic flux, in the first order of Aj_, appears quantized ($*°* = 
or 1 + O(AjJ), for A_|_ < 0.12). This transition is reflected in figure f|. 

An appealing question is that of the behaviour - analysis of the robustness and stability 
specially - of different metastable configurations of a system under the variations of the 
external parameters. It provides information on the energy landscape of a system and 
other properties of the phase space which can determine interesting situations such as a 
constrained dynamics or slow relaxation processes. The inclusion of the CIMF enlarges the 
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stability range of the vortex configurations, as we have explicitly shown in the extreme cases 
of one single vortex and a w = | phase. Quite recently Hwang, Ryu and Stroud0 have 
extended the results on the ground state properties of the Josephson junction ladders in 
the large penetration depth limit, to treat the IV characteristics of a ladder array. They 
found that when "current is injected perpendicular to the ladder edges, the critical current 
is unchanged from its / = value up to a penetration field of f& ~ 0.12 flux quanta per 
plaquette" . This result can be easily understood from the considerations on the stability of 
the one vortex configuration we have made above. At low values of the external field the 
ground state phase configuration of a bidimensional square Josephson junction array has a 
vortex density different from zero and the existence of vortices produces a critical current 
lower than that of the u = phase. In the ladder, however, the situation is quite different: 
at low values of / the ground state in the ladder is the uj = phase and, moreover, as we 
have described above, the configuration with just one vortex is unstable and the u = phase 
is the only stable attractor for arbitrary initial phase configurations. Consequently, at low 
values of the external field the critical current of the ladder is expected not to change, since 
it depends on the vortex density. However, we have shown that the critical value of the 
field for the stability of the one single vortex configuration is = 0.1175 ± 0.00065 for an 
isotropic ladder. For values of the field above different vortex configurations are stable; 
thus, arbitrary initial phase configurations generically relax to different possible metastable 
configurations, which can be described as irregular arrays of vortices. In that situation the 
critical current is essentially associated with the depinning transition of these vortex phases, 
which occurs at a lower value of the external current. Hwang, Ryu and Stroud give a value 
around 0.12, quite close to our prediction of /J = 0.1175 ± 0.00065. 

As we have seen above, as A^ decreases, the value of also decreases, vanishing at 
Aj_ = 1.812 ± 0.018. Then, the diagram of stability (Fig. 0) allows to make the following 
conjecture on the behaviour of the critical current of the Josephson junction ladder when 
current induced magnetic fields are taken into account: as the penetration depth decreases, 
the range of values of the external field for which the critical current remains unchanged0, 
also decreases and it vanishes when A^ = 1.812 ± 0.018. This conjecture is based on the 
natural association between the critical current of the ladder and the stability of the one 
vortex configuration. 
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APPENDIX A: THE INDUCTANCE MATRIX 



The L matrix is obtained in the following form: we have applied the Biot-Savart law in 
order to calculate the magnetic field induced on a link by all the currents circulating in the 
array. We thus obtain 



I d in dl = FF a 0.^sI-yS, (J) 



where FF is the form matrix, which depends on the geometry of the array, and I a p the 
total current passing through link af3. If the links a(3 and 7<5 are perpendicular we take 

FFap-^s = 0. 

The self- inductance term FF a p. a p depends strongly on the form of the current. If, for 
example, it is supposed to flow within a tube of length I with circular cross section of radius 
r it is given by 

FF a0 . af3 = 21 [log - ?J . (8) 

If I is measured in meters, FF is given in 10 _1 Henries. 

Instead, the mutual inductances between different links are sensibly the same as that 
of the filaments through the centers of their cross sections, even when the links are very 
close. In particular, in the case of cylindric currents, the mutual inductances are absolutely 
independent from the radius. 

In order to express |7| as a function of adimensional quantities we introduce // = -^FF, 
thus obtaining 

- ^ f c? nd dl = -L- £ // aft7 ^, (9) 



a af3 



where currents are normalized to the critical current of the link I c (in the case of an 
anisotropic ladder, where I x ^ I y , we take I c as I cx ). X± is the penetration depth, see 
eq. (|). 

It is easy to obtain an equivalent description of the self-field effect in terms of the mesh 
currents and the magnetic flux, as required in equation (^). The magnetic flux on z-cell is 
given by 

*?* = f$ nd dl = p- £ a^ d } (10) 

where the sum is over the four links aj3 of cell i, and can be expressed by a linear equation 
of the form 

£ a% = A, a ^. (11) 
We have used greek and roman symbols to denote, respectively, links and cells. 
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On the other hand, also the mesh currents (zj) are related to the link currents {i a p) 
through a linear operator 

= X] B a f3;iii- (12) 
a/3£i 

Combining (|), (0), (0) and flTJ) we obtain 

$- nd = ^ii^j , ^ij = ~r- Q 2 , , A i;i = X!E^;^«ft7^;j' ( 13 ) 

The Aj.j elements depend on the distance \i — j\ between the cells considered. The general 
properties of matrix are: A o is positive, and A^ < for i ^ j] for two cells far away 
enough (\i — j\ > 10) the mutual inductance is |Ay| ~ \i — j|~ 3 , as calculated ini. 

In the following, we will consider that currents flow inside cylindric tubes with circular 
cross section of radius r. If currents are supposed to have a very small cross section (com- 
pared to the lattice spacing, a), then f f a /3-, a /3 » ffap^s- Then it can be easily verified 
that Ajj ~ 81og(l/r), and Aj i+1 ~ — 21og(l/r). The other terms | A ii+J - 1 (^>i) << A& and 
can be neglected. As the cross section increases, ff a /3;a/3 becomes smaller while the rest 
of the ff's remain sensibly the same; Aj j decreases, A^+i increases and changes sign at 
r ~ 0.24a and A^i + j, (j > 1) remain the same. There is a range of r values r e [0.2,0.28] 
where |A M | > 20\A i:j \(j ^ i). 

Thus we can establish r ranges where the different approximations A,B and C made in 
the paper are acceptable. If r << a, Aj j > |Aj j +1 | >> |Aj j + j|(j > 1) and approximation B 
(considering just the self-inductance plus the nearest neighbours terms in A) is justified. As 
r increases, the terms Ai t i+j(J > 1) also increase, but are yet too small. They give just small 
corrections to the final results. In a narrow range of r values around r ~ 0.25 the dominant 
contribution to the inductance is self-term (and A is a good — th order approximation). 
Finally, for greater r \A iji+ j\/ A iti cannot be neglected and considering the whole inductance 
matrix is compulsory. 

In our calculations, we have considered square cells; for the study of case C we have 
considered an intermediate value of r (currents are supposed to flow within a cylinder of 
length a and ratio 0.005a@). In this case, A 0;0 = 38.194 and ^ = { 1, -0.20332, -0.0040118, 
-0.0010570 ... }. 

We can now calculate the equivalent dx in the case / = | for the configurations u = 
and oj = |. This can be easily made if one considers the spatial periodicity of the solutions. 
If o) = all the cells in the array have the same flux and current; thus we can define an 
effective self-inductance matrix as 



fi = (Lj.j + 2 * (Lj.j+i + Lj.j +2 + ...)) %i = L e ffU, (14) 

where factor 2 is due to the sum of the contributions from the cells on the left and on the 
right. The terms Lij,i ^ j are negative and thus L e jf < La. Now dxeffi^ = 0) = j^— = 
ll,176Aj_. 

In an analogous way, one can calculate the value of L e ff for a u = \ solution: for f = | 
the flux and the current in one cell have the same magnitude and the inverse sign of those 
in the adjacent plaquettes. Thus 
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fi = (Li-i + 2 * (-L i;i+1 + Li-i+2 - ...)) k = L e ffii, (15) 
where now L e ff > La and dxeffi^ = \) = Sti" 3 /-^// = 4, 638Aj_. 
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FIGURES 



FIG. 1. Schematic representation of the Josephson junction array we study here: an anisotropic 
(J x 7^ J y ) ladder, in the presence of an external field. The sites denote superconducting islands 
and the crosses the junctions themselves. Right-most plaquette shows the mesh current ii and the 
gauge choice, here f\ ot = fo + fi where fo = is the flux due to the external field and fi the 
induced flux in the plaquette i. 

FIG. 2. Ground state phase diagrams of the JJL obtained using the method of effective poten- 
tials. Each phase is defined by the value of oj and, for clarity, only a few of the transition lines 
are represented. Figure (a) shows the results for the no screening field case (or A^ — > oo). (b) 
Phase diagram for a A^ = 1.0 ladder using approximation A (diagonal inductance matrix) to the 
calculation of the induced fluxes, (c) Phase diagram for a Aj_ = 1.0 ladder using approximation B 
(self plus nearest neighbours inductances) to the calculation of the induced fluxes. 

FIG. 3. Vortex shape as a function of the penetration depth A^ and the anisotropy. We 
consider an isotropic 128-cell ladder with just one vortex. The whole inductance matrix is used. 
fo = 0. The shape is symmetric, and we just draw the right-half of the vortex (the origin of 
coordinates is the central plaquette of the ladder). The figure shows the gauge- invariant phase 
difference along the horizontal links belonging to the upper branch of the ladder. We compare the 
cases with J y = J x , Aj_ = 1 (black circles) J y = 0AJ x , Aj_ = 1 (squares) and J y = J x , X± = 0.012 
(rhombs). The inset shows the induced flux around the central plaquette, in the same cases as 
before. 

FIG. 4. Induced flux as a function of Aj_ for a configuration u = ^, at fo = \, Jx = Jy We 
draw the flux through two plaquettes with vorticity and 1, respectively. We consider the C case: 
whole inductance matrix. In the inset we show the derivative of the induced flux respect to the 
logarithm of the penetration depth. We can distinguish three regions: the extreme limits A^ > 4 
(the induced flux is fi = + 6/X±) and Aj_ < 0.12 ( fi = rii — fo + 5X±) and an intermediate region 
around Aj_ ~ 0.7, where the derivative is maximum. 

FIG. 5. Energies of different configurations as a function of frustration (e w (fo)) for a penetration 
depth Aj_ = 0.5 in the case A. We study configurations w = 0, 1/5, 1/4, 1/3, 2/5 and 1/2 (in order 
of decreasing slope). We consider an isotropic ladder (J x = J y ). The ground state energy at each 
value of fo is given by the envelope of the curves. The inset shows the value of w corresponding to 
the minimum energy curve for each fo- 

FIG. 6. Devil's staircase observed in case C in the limit Aj_ — ► for an isotropic ladder, 
calculated as explained in the text. 
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FIG. 7. Border between stability and instability regions in the parameter space for the extreme 
cases of one single vortex (black points) and w = 1/2 configuration (rhombs) in an isotropic ladder. 
In the first case we have considered a 128-cell ladder. We have checked that the curves fit to a 
quadratic function: /o = a (j3 + j^jj^yj * ( a]~7o) ~~ \L TToj )" ^ ±c ^ 1S * ne penetration depth 
below which a configuration is stable at /o = 0, and a/3 x ± l(o) = wnere fc is the value of the 
frustration below which a configuration is no more stable in the limit of no inductance (Aj_ — > oo). 
For a configuration with one single vortex, f c = 0.1175 ± 0.00065 and Aj_ c (0) = 1.812 ± 0.018; in 
the w = 1/2 case, f c = 0.215 ± 0.001 and A ±c (0) = 1.197 ± 0.006. 
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